Investigation of the monopole magneto-chemical potential in spin ices using capacitive torque magnetometry

The single-ion anisotropy and magnetic interactions in spin-ice systems give rise to unusual non-collinear spin textures, such as Pauling states and magnetic monopoles. The effective spin correlation strength (Jeff) determines the relative energies of the different spin-ice states. With this work, we display the capability of capacitive torque magnetometry in characterizing the magneto-chemical potential associated with monopole formation. We build a magnetic phase diagram of Ho2Ti2O7, and show that the magneto-chemical potential depends on the spin sublattice (α or β), i.e., the Pauling state, involved in the transition. Monte Carlo simulations using the dipolar-spin-ice Hamiltonian support our findings of a sublattice-dependent magneto-chemical potential, but the model underestimates the Jeff for the β-sublattice. Additional simulations, including next-nearest neighbor interactions (J2), show that long-range exchange terms in the Hamiltonian are needed to describe the measurements. This demonstrates that torque magnetometry provides a sensitive test for Jeff and the spin-spin interactions that contribute to it.

1 The (0 0 1) rotational plane Supplementary Figure 1 shows the schematic of the capacitive torque magnetometry (CTM) set-up and a 3D representation of the rotation planes. It also highlights the magnetic phases that we experimentally observe in both rotational planes when the magnetic field is aligned with the various high symmetry axes that lie on these planes. Based on the rotation schematics shown in Supplementary Figure 1 with B the strength of the applied field. When rotating in the (0 0 1) plane, the zero crossings are observed at intervals of 45 • and sharp turnovers in the data occur near any of the ⟨1 1 0⟩ family of directions. The system always remains in one of the Pauling states in this plane, in particular (2:2) 0 (or Q = 0 1 ), with the field direction determining which of the 6 equivalent Pauling states dominates. These states are characterized by a saturated moment pointing along either thex orŷ orẑ directions. For our chosen plane of rotation, it is either alongx orŷ. We observe a sharp turnover in the data with a zero torque crossover points at or near the ⟨1 1 0⟩ directions. This is when all β -spins flip and m sharply rotates by 90 • . The intermediate phases that occur near the ⟨1 1 0⟩ directions could be understood as an equal volume fraction mixture of (2:2) 0 phases with half of the sample having m ∥x and other half of the sample having m ∥ŷ. This would result into a measured saturation magnetization of about 4.1 µ B / Ho 3+ which is same as for the (2:2) X phase. However, in this rotation plane, we do not see any evidence of a stable (2:2) X phase as that would require a long range ordering among the β -spins. Table 1. Functional forms for calculated torque curves for the magnetic spin textures associated with the (0 0 1) rotational plane that are shown as solid/dotted curves in Figure 3 of the main text. Note, the magnetization, m, is given as the magnetization per one unit cell, i.e., for 16 Ho 3+ sites, with µ=10µ B per site. Transitions between the indicated spin textures requires reversal of 4 β -spins, as indicated in red. For each of the mentioned phases the torque response in a given applied field as a function of angle θ , defined as the angle away from the [  the associated net magnetic moments, and the functional forms of the torque curves as B rotates in the (0 0 1) plane. Transitions between spin textures, which can involve spin flips on either the α or β sublattice, lead to gradual changes in the value for m = µ √ 3 · S, with µ = 10µ B the moment per Ho 3+ ion and with S = (S x , S y , S z ). In the (0 0 1) plane, transitions involve only β -spin flips and each flip leads to ∆S x = ±2 and ∆S y = ±2 while S z remains zero. The torque curves for these intermediate spin textures are calculated in the same way as the other curves, taking into account only this change in S (see Figure 3 of the main text for the calculated torque response as a function of B(θ )). The field-angle phase diagram for the (0 0 1) rotational plane as determined from the CTM measurements, is provided in Supplementary Figure 2  In the plane that contains the body diagonal [1 1 1], a similar procedure is used to calculate the angular dependence of the torque response. For this plane, shown in Supplementary Figure 1 (d) and with B ≥ 2 T, the torque response shows additional zero crossings and turnovers with significant changes in the observed torque amplitude when approaching different high symmetry axes. Note, in this rotational plane α-spins are involved in the phase transition between (2:2) 0 and the (  and β -spins are involved in the phase transition between the (2:2) X and (3:1) monopole phases. We define Table 2. Functional forms for torque curves for the magnetic spin textures associated with the (1 1 0) rotational plane that are shown as solid curves in Figure 3 of the main text. Note, the magnetization, m, is given as the magnetization per one unit cell, i.e., for 16 Ho 3+ sites, with µ=10µ B per site. Transitions between the indicated spin textures requires reversal of 4 α-or β -spins, as indicated in blue and red text, respectively.
Phase    leads to ∆S z = ±2 while S x and S y remain unchanged. After every α-or β -spin flip in a unit cell consisting of four tetrahedra, the system goes into a mixed domain phase where 25% of the system has made the phase transition (i.e., one tetrahedron in every unit cell has transitioned between a Pauling and a monopole state). Therefore, the Pauling:monopole volume fraction changes from 4:0 → 3:1 → 2:2 → 1:3 → 0:4 during a phase transition. This provides the appropriate scaling factors to extract the magnetic moment m, and hence, the functional form of the torque equations for the intermediate phases during the phase transition. The torque curves for the ordered phases and these intermediate spin textures are shown in Figure 3 in the main text as solid and dotted curves, respectively. The field-angle phase diagram for the (1 1 0) rotational plane as determined from the CTM measurements, is provided in Supplementary Figure 2

Integrated torque curve and the magnetic anisotropy energy
From the measured torque response we can extract the anisotropy energy (E a ) as a function of field strength and direction. The magnetic torque is the derivative of this energy, where θ is the rotation angle around a defined axis within the crystal. Measuring torque responses for different rotation planes thus provide insight into the shape of the anisotropy energy surface. In most systems this torque provides a direct means to determine the magneto-crystalline anisotropy, as most magnetic materials have a rather small crystal field that can be easily overcome by an applied field. Since, the spin-ice system has a very large crystal field splitting 2 , rather than measuring the magneto-crystalline anisotropy surface, we are probing the energy surface associated with the various noncollinear spin textures that are stabilized under different magnetic field conditions. Integrating the torque curves reveal the relative energy differences between these possible spin textures as the sample is rotated in applied fields. Field-scaled integrated torque curves for both the (0 0 1) and (1 1 0) planes are plotted in Supplementary Figure 3 for various applied field strengths. We assume that the energy surface is continuous. Thus, it becomes clear that the (2:2) X phase, which is expected to form when the field is applied along any of the ⟨1 1 0⟩ directions, resides on a sharp energy maximum (i.e., is not stable) when rotating in the (0 0 1) plane, yet it resides on a local minimum in the energy surface when rotating in the (1 1 0) plane. It appears to reside on a saddle point in energy. Based on this energy surface, a long-range ordered (2:2) X phase does not appear to form when the field is confined to rotate in the (0 0 1) plane, yet an ordered (2:2) X phase does form when the field is confined to rotate in the (1 1 0) plane, with the angular range of this phase clearly showing a strong dependence on field-strength. direction. We express the applied field as B = B[cos(θ )ĵ ′ + sin(θ )k ′ ] and τ n =n · τ = [1,1,0] √ 2 · (m × B), identical to the reference frame used for the (1 1 0) rotation. We note that θ , defined as the angle between B and the [1 1 2] direction, is fixed during the measurement (θ = 90 • -∆θ ), while B varies between 7 T and -7 T. Supplementary Table 3 lists all the moment vectors and the torque curves for the stable spin textures. It also indicates which and how many spins per unit cell are flipped during a phase transition. When a β -spin flips, each flip leads to ∆S z = ±2, while S x and S y remain unchanged whereas, an α-spin flip leads to ∆S x =∆S y = ±2 with S z remaining unchanged.
Assuming ∆θ = 5 • , the experimentally measured ratio of the slopes between the (3:1) phase and the (2:2) X phases gives a value of 18. This is extremely close to the expected value of 22 based on magnetization measurements reported by others 3,4 . An agreement within 5 % between the expected and experimental values gives us a great deal of confidence in the measurement set-up and phenomenological model. We also recognize that there is no stable phase when the field reverses its orientation. This occurs between -0.1 T → -1.8 T and 0.1 T → 1.8 T when all α-spins flip. However, approaching the same region from the high field end (±7 T) results into a clean (2:2) X phase, which follows after the β -spins flip. Because different sets of spin-flips are involved during the low and high field transitions in the field sweep measurement, a hysteretic butterfly loop within the field range of -1.8 T< B <1.8 T results. The energy cost for a spin flip to occur and create the (3:1) state on a single tetrahedron is likely too high at low field, i.e., when the (3:1) state is inaccessible, hysteresis appears, with the 8-α spins likely flipping in pairs. Table 3. Functional forms for torque curves for the magnetic spin textures associated with the field sweep measurements. For the solid curves in Figure 4(b) of the main text, with θ = 85 • . The magnetization, m, is given as the magnetization per one unit cell, i.e., for 16 Ho 3+ sites, with µ=10µ B per site. Transitions between the indicated spin textures requires reversal of 8 α-or 4 β -spins, as indicated in blue and red text, respectively.  Figure 3 (a) in the main text, the steep transition around the ⟨1 1 0⟩ family of directions develops obvious hysteresis at low applied fields at T = 500 mK. The hysteresis is a clear sign of the glassiness of the system. The field rotation speed may be important in determining the degree of hysteresis, as the system will be effectively out of equilibrium if the field rotates too fast. Here we note that hysteresis is also present at higher fields, but gradually closes with increasing field. We have also investigated how the hysteresis changes with increasing temperature (for B = 2 T), Supplementary Figure 4(a) shows the torque response at both T = 500 mK as well as T = 1.6 K. The wine and blue solid lines correspond to calculated torque curves associated with the (2:2) 0 and (2:2) X phases, respectively (see Supplementary Table 1). It is clear that the hysteresis disappears with increasing temperature. Therefore, we associate the appearance of the hysteresis with the local internal fields (spin-spin interactions) of the spin-ice state as it appears below the spin freezing temperature.

The (1 1 0) rotational plane
As mentioned in the main text and as shown in Supplementary Figure 4(b), the transitions around the [1 1 2] and [1 1 2] directions between the (2:2) 0 and (3:1) monopole phase disappear below B = 2 T, i.e., the (3:1) monopole phase does not form at low field. As a result hysteresis appears around the [1 1 1] and [1 1 1] directions for rotations in B = 1 T. This hysteresis is clearly associated with the spin-ice phase as it has fully disappeared at T = 1.7 K (the solid magenta curve in Supplementary Figure 4(b)). The wine, blue, and green solid lines correspond to calculated torque curves associated with the (2:2) 0 , (2:2) X and (3:1) monopole phases, respectively (see Supplementary Table 2). We associate the appearance of the hysteresis with the local internal fields (spin-spin interactions) of the spin-ice state and assume that in this rotational plane it also appears below the spin freezing temperature.
In Supplementary Figure 4(c) we show the CTM response measured at 1.7 K in various applied fields. Using the same procedure as described in the main text we compare these angular sweeps with the calculated torque curves and extract the critical angles for the (2:2) 0 ⇔ (3:1) and (2:2) X ⇔ (3:1) transitions (see Supplementary Figure 4(d)). The critical angles as a function of applied field are plotted in Supplementary Figure 4(e). The thermal smearing of the transitions and the fact that at 2 T we do not observe the (3:1) state, makes the analysis as described in the main text difficult for the high temperature data.

6/12
Instead of fitting the critical angles, we plot them together with the functional forms for B(θ ) using 1.4 and 2.2 K for J α e f f and J β e f f , respectively. The agreement between B(θ ) and the data is very good.

Temperature dependent hysteresis in field sweep with B near the [1 1 1] direction
We have also investigated how the field sweep changes as the temperature is gradually increased from T = 500 mK to T = 1.6 K (see Supplementary Figure 4(f)). For these measurements, the field is misaligned away from the [1 1 1] direction by a small angle of ∆θ ≈ 5 • towards the [1 1 0] direction. While the linear responses at high fields (B > 2 T) and at low fields (0 T < B < 1.5 T) remain, the sharp turnover around B = 1.8 T gradually smears out at higher temperature. The hysteresis around zero field disappears fully when the sample is heated to above T = 0.7 K, i.e., around the spin freezing temperature. We observe no hysteresis around the sharp feature at B = 1.8 T, which is expected and consistent with reports by others 3 . We also notice that as the temperature is increased, the value for ∆C at the turnover (B = 1.8 T) decreases in amplitude. This is consistent with the formation of thermal defects in the (2:2) X phase (with S = (8, 8, 0)) resulting in a gradual increase in S z as the long-range antiferromagnetic order on the β -spins gets randomized. Additionally, there is a drift in the measurement at 1.6 K (a little drift is observed at 1 K as well) and we associate this trend with a systematic drift in the sample temperature during field sweep. So, there is a clear evidence that the capacitance picks up an additional term because of the temperature variation during the field sweep sequence. This is expected in a sorption-pump based He3 cooling system because higher field burns helium at a faster rate.

Hamiltonian and Monte Carlo Simulations
The Hamiltonian of interest for Ho 2 Ti 2 O 7 is, whereS i are classical spin vectors with |S i | = 1. The tilde is used to indicate that the spins are constrained to point along the local ⟨1 1 1⟩ axis of the tetrahedra they belong to. r i is the real-space location of site i, r i j ≡ r i − r j , ⟨i, j⟩ (⟨⟨i, j⟩⟩) refers to nearest-neighbor (next nearest-neighbor) bonds, r nn is the nearest-neighbor bond distance, J 1 (J 2 ) is the nearest-neighbor (next-nearest neighbor) interaction strength, and D is the strength of the long-range dipolar term. The i > j notation guarantees each of pair of spins is only counted once. gµ B is the size of the magnetic moment and B is the applied magnetic field. The spins are effectively Ising-like since they are constrained to be aligned along a local ⟨1 1 1⟩ axis, i.e., where σ i = ±1 is a scalar and z a(i) is a unit vector depending on the sublattice a = 0, 1, 2, 3 to which the site i belongs. The directions of the z a vectors are shown in Supplementary Table 4.
When the dipolar interaction is truncated to the nearest neighbor term, the Hamiltonian becomes, where J s−DSI 1,e f f ≡ J 1 +5D

3
. For the case of zero applied field, J s−DSI 1,e f f > 0 results in the "2 in-2 out ice rule" for the ground state manifold of the spin ice.
We perform classical Monte Carlo (MC) simulations for both the nearest-neighbor (NN) and the standard dipolar spin-ice (s-DSI) models, with the parameters shown in Supplementary Table 5. The s-DSI model includes long-range dipolar interactions via the Ewald summation method [5][6][7][8] , the parameter set for the s-DSI model is taken from Refs. 9,10  1,e f f are provided as a reference, but they are not used in the simulations.
Ref 11 3.41 -0.14 0.025 (0.025)  In Supplementary Table 6 we compare our values to others reported for HTO and Dy 2 Ti 2 O 7 (DTO). We simulate finite size pyrochlore clusters with periodic boundary conditions with N spins = 16 × 4 3 = 1024 lattice sites and employ the Ewald summation method [5][6][7][8] . Demagnetization effects 6 are found to be significant especially for field directions close to the observed phase boundaries. We model this with the help of a surface term in the Ewald summation 7 , Incorporating this term with µ ′ = 1 provides a good match to the experimental observations. Our MC sampling procedure employs the standard accept-reject Metropolis algorithm with an equal probability of "loop" moves (effective for sampling the low energy spin-ice configurations at low temperature) 6,16,17 , and single spin flip moves (effective at high temperature). We measure the thermal average of the magnetization ⟨m⟩ (per unit cell) using 10 5 MC samples. One sample per MC sweep is recorded for thermal averaging, and one MC sweep comprises of a set of N spins proposed moves (which are individually of loop or single spin flip type). We then compute the component of the torque alongn (per unit cell) τ n ≡n · (⟨m⟩ × B), as discussed in the previous section. This quantity corresponds directly to what is being measured in the torque magnetometry experiment (up to device specific constants). We scan the space of field alignment angles in steps of 1 • (or smaller near phase boundaries), considering the same range of angles as that recorded in the experiment.
In Supplementary Figure 5 we show snapshots of spin textures recorded during the MC simulations at T = 0.5 K as the field was rotated from the [001] direction at 0 • . We show snapshots at various relevant angles, observing the expected stable phases (2:2) 0 at 0 • , (3:1) at 45 • and 100 • , and (2:2) X at 50 • and 51 • . The reason that two (2:2) X configurations are shown is because in the 47 -63 • range, the spin texture toggles between the two possible states with oppositely oriented β -spin chains every couple of degrees, indicating that these textures are indeed degenerate. No defects are observed in either of the (2:2) X configurations, suggesting that thermal fluctuations at T = 0.5 K are small and the state is relatively stable. Lastly, the configurations taken at angles of 4 • , 46 • , and 64 • show intermediate or transient spin textures in which the two stable phases involved in the transition form a domain structure. This is visible from the presence of the double arrows, which indicates that in that column of spins, both in and out spins (either unflipped α-or β -spins) exist.
In Supplementary Figure 6, simulated torque curves are compared to the experimental curves taken at 0.5 K in 2 and 11 T, respectively. The results for a strictly nearest neighbor model (NN, blue curve) and for the s-DSI (long-range dipolar) model (including Ewald summation, loop moves, and demagnetization effects 6 , red curve) are shown. The wine, orange, green, and purple curves were simulated with a more generalized DSI model using s-DSI model parameters and including various values for the J 2 interaction (0.01, 0.02, 0.03, and 0.04 meV).
where ⟨⟨⟨i, j⟩⟩⟩, a and ⟨⟨⟨i, j⟩⟩⟩, b refer to third nearest neighbors of a and b types respectively. (The a type third nearest neighbors share a common neighbor and are thus connected to each other via a minimum of two "hops" on the pyrochlore lattice. The b type third nearest neighbors do not share a nearest neighbor, they are instead related by a path involving a minimum of three hops on the pyrochlore lattice 11,15 .) We enumerated the first, second and third nearest neighbors (of both kinds a, b) for the pyrochlore lattice. (This was done numerically for a cubic unit cell of 16 × 2 3 = 128 sites). We computed the Zeeman energy and the interaction energy per 16-site unit cell for three phases: (2:2) 0 (with m||z), (2:2) X , and 3:1 (+ − −−, i.e., 3-in 1-out).
The sum of the interaction and Zeeman contributions (i.e., the second term in Supplementary Equation 6 with B defined as (h x ,h y ,h z )) gives us the total energy (which per 16-site cubic unit cell) is as below: The magnetic field components as a function of the angle (θ ) are (see Supplementary Equations 2(a) and (b)): To obtain the phase boundary of the (2 : 2) 0 ⇔ (3 : 1) transition at zero temperature (T = 0) we equate the total energies of the two phases (per 16-site unit cell). We get, The phase transitions between symmetry related phases (2:2) 0 (with m|| − z) and 3:1 (+ + +−, i.e., 3-out 1-in) can be evaluated in a similar way. Note that the J a 3,e f f and J b 3,e f f terms cancel out for this transition. Thus, J 3 interactions do not affect this transition.
The effective parameters can be expressed in terms of the original generalized DSI model (now incorporating second and third nearest neighbors with couplings J 2 , J 3a , and J 3b ) following Ref. 15 , who considered the case of J 2 = J 3a = J 3b = 0, but with the opposite sign convention of J 1 compared to our work.The relations are, The factor of 1/3 arises from the geometric factor z i · z j = −1/3 for two spins on different sublattices as is the case for nearest and second nearest neighbors. The third nearest neighbors are on the same sublattice for which z i · z i = 1. Note that this model truncates the dipolar interaction, and is thus expected to be only qualitatively accurate. Based on Supplementary Equations 11 and 10 we can extract a J α e f f that would correspond to this combination of interaction terms.
J α e f f ≡ J 1,e f f + 2J 2,e f f = For J 1 = −0.13499 meV (-1.56 K), D = 0.121567 meV (1.41 K), and J 2 = 0.030 meV (0.35 K) we get J α e f f = 1.52 K, which is lower than J 1,e f f =1.83 K, encouragingly it is in the same ballpark as that measured in the CTM experiment. Hence, once the interaction terms are known, we can rearrange Supplementary Equation 10 and plot B as a function of θ , which represents the phase boundary at which the transition occurs (see Supplementary Figure 7 (a)). We plot the phase boundary for three combinations of J 1 , J 2 , and D.